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Abstract 

We study the classical XY (plane rotator) model at the Kosterlitz- 
Thouless phase transition. We simulate the model using the single 
cluster algorithm on square lattices of a linear size up to L = 2048. 
We derive the finite size behaviour of the second moment correlation 
length over the lattice size £,2nd/L at the transition temperature. This 
new prediction and the analogous one for the helicity modulus T are 
confronted with our Monte Carlo data. This way (5kt = 1.1199 is 
confirmed as inverse transition temperature. Finally we address the 
puzzle of logarithmic corrections of the magnetic susceptibility x & t 
the transition temperature. 



PACS numbers: 75.10.Hk, 05.10.Ln, 68.35.Rh 



1 Introduction 



We study the classical XY model on the square lattice. It is characterised 
by the action 

S = -I3^s x s x+P , , (1) 

X,(l 

where s x is a unit vector with two real components, x = (xi,X2) labels the 
sites on the square lattice, where X\ G {1, 2, L{\ and x 2 £ {1, 2, L 2 } *, 
\i gives the direction on the lattice and jl is a unit-vector in the /z-direction. 
We consider periodic boundary conditions in both directions. The coupling 
constant has been set to J = 1 and (5 is the inverse temperature. In our 
notation, the Boltzmann-factor is given by exp(— S). Sometimes in the liter- 
ature the present model is also called "plane rotator model" , while the name 
XY-model is used for a model with three spin-components. 

Kosterlitz and Thouless [1] have argued that the XY-model undergoes a 
phase transition of infinite order. The low temperature phase is characterised 
by a vanishing order parameter and an infinite correlation length £, associated 
with a line of Gaussian fixed points. At a sufficiently high temperature 
pairs of vortices unbind and start to disorder the system resulting in a finite 
correlation length £. In the neighbourhood of the transition temperature 
T KT it behaves as 

C-aexp^r 1 / 2 ) , (2) 

where t — (T — Tkt)/Tkt is the reduced temperature and a and b are non- 
universal constants. In subsequent work (e.g. refs. [2, 3]) the results of 
Kosterlitz and Thouless had been confirmed and the arguments had been 
put on a more rigorous basis. 

This rather good theoretical understanding of the Kosterlitz-Thouless 
(KT) phase transition is contrasted by the fact that the verification of the 
theoretical predictions in Monte Carlo simulations had often been inconclu- 
sive or even in contradiction. Only starting from the early nineties, Monte 
Carlo simulations allowed to favour clearly the KT-behaviour (2) over a power 
law £ oc t^", which is characteristic for a second order phase transition. A 
typical example for such a work is ref. [4], where the XY model with the 
Villian action [5] was studied on lattices of a size up to 1200 2 . 

*In our simulations we use L\ = L 2 = L throughout 
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The difficulties in Monte Carlo simulations might be explained by loga- 
rithmic corrections that are predicted to be present in the neighbourhood of 
the transition. 

In the present paper we like to address two puzzling results presented in 
the literature that are related to this problem: 

• The two most precise results [6, 7] for the transition temperature Tkt 
of the XY-model differ by about 8 times the quoted errors. 

• The magnetic susceptibility is predicted to scale as x °c L 2 ~ n '(In L)~ 2r 
with 7] — 1/4 and r = —1/16 at the transition temperature, t However 
the authors of refs. [10, 11] find in their Monte Carlo simulations 
r = -0.023(10) * and r = -0.0270(10), respectively. 

In refs. [13, 7] the authors have shown that XY models with different 
actions share the universality class of the BCSOS model. This had been 
achieved by matching the renormalization group (RG) flow of the BCSOS 
model at the critical point with that of the exact duals [14] of the XY models 
using a particular Monte Carlo renormalization group method. As a result 
of this matching the estimate f3 KT = 1.1199(1) = 1/0.89294(8) for the XY 
model (1) has been obtained. § The BCSOS model is equivalent with the 
six- vertex model [15]. The exact result for the correlation length of the six- 
vertex model [16, 17, 18] shows the behaviour of eq. (2) predicted by the KT- 
theory. The main advantage of the matching approach is that the logarithmic 
corrections and in particular also subleading logarithmic corrections are the 
same in the XY-model and the BCSOS model. * 

In a more standard approach, Olsson [6] and Schultka and Manousakis 
[19] have studied the finite size behaviour of the helicity modulus arriving at 
the estimates l/f3 KT = 0.89213(10) and l/(3 KT = 0.89220(13), respectively. 

^Note that the analogous result \ oc £ 2 ~' ) (ln£;)~ 2r for the thermodynamic limit in the 
high temperature phase does not hold. In refs. [8, 9] it was argued and numerically verified 
that instead x K £ 2 ? '(1 + c /( m £ + u ) 2 + •••) is correct. 

*The authors confirmed their numerical result for r by a study of Lee- Yang zeros [12] 
§In the case of the Villian action, the matching method gives Pv.kt = 0.7515(2), while 
the authors of ref. [4] had found Pv,kt = 0.752(5) fitting their data for the correlation 
length with the ansatz (2) and a similar fit for the magnetic susceptibility, 
brief discussion of this fact will be given in section 3. 
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These authors studied lattice sizes up to L = 256 and L = 400, respectively. 
While in their approach leading logarithmic corrections are taken properly 
into account, sub leading logarithmic corrections are missed. This might ex- 
plain the missmatch of the results for the transition temperature. Here we 
shall resolve this discrepancy by brute force: We study the helicity modulus 
(and in addition the second moment correlation length) on lattices up to 
L = 2048. 

Having an accurate estimate of Tkt and numerical results for large lattice 
sizes at hand, we then study the scaling of the magnetic susceptibility. Here 
it turns out that the puzzling result for the value of the exponent r can be 
resolved by taking into account subleading corrections. 

A major purpose of the present paper is to check the reliability of standard 
methods to determine the temperature of the transition and to verify its 
KT-nature. This aims mainly at more complicated models, e.g. quantum 
models or thin films of three dimensional systems with nontrivial boundary 
conditions, where the duality transformation is not possible, and hence the 
method of refs. [13, 7] can not be applied. 

The outline of the paper is the following: In the next section we give 
the definitions of the observables that are studied in this paper: the helicity 
modulus, the second moment correlation length and the magnetic suscepti- 
bility. Next we summarise some results from the literature on duality and 
the RG-flow at the KT-transition. We re-derive the finite size behaviour of 
the helicity modulus at the transition temperature. Along the same lines we 
then derive a new result for the dimensionless ratio ^nd/L. This is followed 
by Monte Carlo simulations using the single cluster algorithm for lattices of 
a linear size up to L = 2048 for /3 = 1.1199 and /3 = 1.12091. Fitting the 
data for (3 — 1.1199 we find the behaviour of the helicity modulus and ^2nd/ L 
predicted by the theory for the transition temperature, while for (3 = 1.12091 
there is clear missmatch. Finally we analyse the data of the magnetic sus- 
ceptibility at = 1.1199. 
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2 The observables 

In this section we shall summarise the definitions of the observables that we 
have measured in our simulations. The total magnetisation is defined by 

M = £4 • (3) 

X 

The magnetic susceptibility is then given as 

X = ^M 2 • (4) 

2.1 The second moment correlation length ^2nd 

The second moment correlation length on a lattice of the size L 2 is defined 
by 

p 1 fx X /2 ~ 

Ud ~ 2sin(7r/L) [f " \) ' (5) 
where x is the magnetic susceptibility as defined above and 

1 

L2 



F = — Yji^xSy) cos(27r(yi - x x )/L) . (6) 



Note that the results obtained in this paper only hold for the definition of 
£,2nd given in this subsection. 

2.2 The helicity modulus T 

The helicity modulus T gives the reaction of the system under a torsion [20]. 
To define the helicity modulus we consider a system, where rotated boundary 
conditions in one direction are introduced: For pairs x, y of nearest neighbour 
sites on the lattice with x± — L 1; y 1 — 1 and x 2 = y 2 the term s x s y is replaced 
by 

s x -R a s y = 4 1 ) (cos(«)4 1 ) + sin(a)4 2 ))+4 2 ) (cos(a)4 2 ) - siii^)^) ■ (7) 

The helicity modulus is then defined by the second derivative of the free 
energy with respect to a at a = 

Li d 2 \nZ(a) 

a=0 



T = 

L 2 da 2 
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Note that we have skipped a factor one over temperature in our definition 
of the helicity modulus to obtain a dimensionless quantity. It is easy to 
write the helicity modulus as an observable of the system at a = [21]. For 
Li = L 2 = L we get 

3 KT-theory 

In this section we summarise results from the literature that are relevant for 
our numerical study and also derive a novel result for the finite size behaviour 
of the second moment correlation length at the transition temperature. 

XY models can be exactly mapped by a so called duality transformation 
[14] into solid on solid (SOS) models. E.g. the XY model with the action (1) 
becomes 

{h} x,n 

where the I n are modified Bessel functions and the h x are integer. The XY 
model with Villian action [5] takes a simpler form under duality: 

^°* = £exp(-^E(^-W 2 ) , 

{h} \ Z ' J x,n ) 

where the h x are integer again. This model is also called discrete Gaussian 
(DG) model. In the context of finite size scaling one should pay attention to 
the fact that the boundary conditions transform non-trivially under duality. 
E.g. periodic boundary conditions in the XY model require that in the SOS 
model one sums over all integer shifts hi and hi at the boundaries in 1- and 
2-direction, respectively. 

It turned out to be most convenient to study the Kosterlitz-Thouless 
phase transition using generalisations of SOS models (see e.g. refs. [2, 3]). 

3.1 The Sine-Gordon model 

The Sine-Gordon model is defined by the action 

S SG = -^E COS ( 27T ^) > ( 12 ) 
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where the variables 4> x are real numbers. For positive values of z, the periodic 
potential favours <p x close to integers. In particular, in the limit z — > oo, we 
recover the DG-SOS model. In the limit z = we get the Gaussian model 
(or in the language of high energy physics, a free field theory). The Sine- 
Gordon model (using cutoff schemes different from the lattice) can be used 
to derive the RG-flow associated with the KT phase transition. For (3 > 2/tt 
the coupling z is irrelevant, while for (3 < 2/ir it becomes relevant. To discuss 
the RG-flow it is convenient to define 



The flow-equations are derived in the neighbourhood of (x,z) = (0,0). To 
leading order they are given by 



where t — In I is the logarithm of the length scale / at which the coupling is 
taken. Note that we consider a fixed lattice spacing and a running length 
scale /, while e.g. in ref. [3] the cutoff scale is varied. This explains the 
opposite sign in the flow equations compared with e.g. ref. [3]. The const 
in the equation above depends on the particular type of cut-off that is used. 
Corrections of 0(z 3 ) have been computed in ref. [3] and confirmed in ref. [22]. 
Here we are mainly interested in the finite size behaviour at the transition 
temperature. Therefore the trajectory at the transition temperature is of 
particular interest. It is characterised by the fact that it ends in (x, z) = 
(0, 0). To leading order it is given by 



x = tt(3 - 2 . 



(13) 



dz 

di 

dx 

di 



xz + 



const z 2 + ... , 



(14) 



(15) 



x = const 1 ! 2 z _ 



(16) 




(17) 



1 



(18) 



x = 



\nl + C 
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where C is an integration constant that depends on the initial value Xi of x 
at I — 1. Taking into account the next to leading order result of ref. [3] the 
flow on the critical trajectory becomes 

1=-*"- ¥ - ■ < 19 > 

Implicitly the solution is given by [3] 

, , 1 1 1, l/x + 1/2 

In = ln^ '— , 20 

x Xi 2 l/xi + 1/2 ' V ; 

where now the initial value x^ of x takes the role of the integration constant. 
The authors of ref. [3] give an approximate solution of this equation that is 
valid for Xi » x. This leads to corrections to eq. (18) that are proportional 
to In | lnL|/| lnL| 2 . However, in our numerical simulations we are rather in 
a situation where Xi and x differ only by a small factor. Therefore we make 
no attempt to fit our data taking explicitly into account the last term of 
eq. (20). 

An important result of ref. [3] is that corrections proportional to 
In | lnL|/| lnL| 2 arise from the RG-flow in the (x, z)-plane and are not caused 
by some additional marginal operators, which might have different ampli- 
tudes in different models. Therefore the two-parameter matching of refs. 
[7, 13] is sufficient to take properly into account corrections proportional to 
In | lnL|/| lnL| 2 (and beyond). 



3.2 Finite size scaling of dimensionless quantities 

Here we compute the values of the helicity modulus T and the ratio (,2nd/ L 
at T KT in the limit L — > oo and leading 1/lnL corrections to it. Since for 
both quantities the coefficient of the order z is vanishing, this can be achieved 
by computing both quantities at z = (i.e. for the Gaussian model) and 
plugging in the value of (5 given by eq. (18). 



3.2.1 The helicity modulus 

The helicity modulus can be easily expressed in terms of the SOS model dual 
to the XY model: 

T = ^(hDsos , (21) 
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where hi is the shift at the boundary in the 1-direction. In this form we 
can compute the helicity modulus in the Sine-Gordon model. To this end we 
have to compute the free energy as a function of the boundary shifts hi,h 2 : 



F(h 1 ,h 2 ) — —\n(Z(h 1 ,h 2 )/Z(0,0)) , 



(22) 



where Z(hi,h 2 ) is the partition function of the system with a shift by hi 
and h 2 at the boundaries in 1 and 2-direction, respectively. From the SG- 
action (12) we directly read off that F(h\, h 2 ) is an even function of z. Hence 
the leading z dependent contribution is 0(z 2 ). Hence for our purpose the 
purely Gaussian result z = is sufficient. For the action (12) at z = we 
get 

Z(h u h 2 ) = fn[<j>\ exp^-^^^-^+A"^) 2 ) 

= |D[0] expf-iL 1 L 2 K + ^) + £(0 ;c _0 



exp {-—^L^dj + dj) 



exp \~W 



Li\ Li 2 




(23) 



where we have defined d^ = h^/L^. Note that we have distributed the 
boundary shift along the lattice by a reparametrisation of the field: 



4>x = 4>x + xidi + x 2 d 2 , 
where 4> x is the original field. It follows 

L 2 E ftl exp(-^^) h\ 



T 



L ^ £ fcl exp(-ig/i?) 



(24) 



(25) 



Alternatively we might evaluate the helicity modulus in the spin-wave 
limit of the XY model on the original lattice. This is justified by the duality 
transformation presented in ref. [2] in appendix D. Here we are only inter- 
ested in the Gaussian limit of the model. Under duality the f3 of the Gaussian 
model transforms as $ — 1/(3. Secondly we have to take into account that 
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even though vortices are not present in the limit z — 0, the periodicity of the 
XY model has to be taken into account for the boundary conditions. Hence, 
the proper spin-wave (SW) description of the XY-model on a finite lattice 
with periodic boundary conditions is 



Zsw= E W(n 1 ,n 2 )Z(0,0) , 



(26) 



m,n 2 



where n\ and n 2 count the windings of the XY-field along the 1 and 2 direction 
respectively. In the Gaussian model they are given by shifts by 27rni and 27rn2 
at the boundaries. The corresponding weights are 



W(ni, n 2 ) = exp 



2/3 



J 2 2 



nf + 



-Tin 



(27) 



Here we can easily introduce a rotation by the angle a at the boundary: 

/ (27T) 2 



Z S W,a = eX P 



ni,n 2 



2(3 



L2 <n 1 + a/(2n)f+ Ll ^ 



Z(0,0) . (28) 



_L 1 L 2 
Plugging this result into the definition (8) of the helicity modulus we get 



T 



L 2 Eni eX P ( 


(271711 ) 2 L 2 \ 


2im\ L 2 


2 


2/3 Lj 







P Ll E„, exp 



(- 



(2ttwi) 2 L 2 
2/3 Li 



) 



(29) 



In the literature often only T — 1/(3 — (3 is quoted and the (tiny) correction 
due to winding fields is ignored. We have checked numerically that the results 
of eq. (25) and eq. (29) indeed coincide. Here we are interested in the case 
of an L 2 lattice in the neighbourhood of (3 — 2/ir. One gets 



T L 2 iZ=0 = 0.63650817819... + 1.001852182... {(3 - 2/tt) + .. 



(30) 



Plugging in the result (18) and identifying the lattice size L with the scale 
at which the coupling is taken, we get 



T L 2 , transi tion = 0.63650817819... + 



0.318899454... 
InL + C 



+ ... 



(31) 



Contributions of 0(z 2 ) that we have ignored here are proportional to l/(mL+ 
C) 2 at the transition. 
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3.2.2 The second moment correlation length 

In this section we derive a result for the dimensionless ratio ^ 2n d/ L analogous 
to eq. (31) for the helicity modulus. To this end we have to compute the XY 
two-point correlation function as a series in z. For the limit L — > oo, the 
result can be found in the literature. It is important to notice that similar to 
the helicity modulus 0(z) contributions to the correlation function vanish. 
I.e. also here the Gaussian result is sufficient for our purpose. The non- 
trivial task is to take properly into account the effects of periodic boundary 
conditions on the finite lattice. The starting point of our calculation is the 
spin wave model (26). Following the definition (24), a difference of variables 
<f> x and (fiy of the system with shifted boundary conditions can be rewritten 
in terms of the system without shift: 

(fix - 4> v = <\>x ~ <\> y + Pin 1 (x 1 - yi) + p 2 n 2 (x 2 - y 2 ) (32) 

with pi = 2n/ Li. Using this results, the spin-spin product can be written as 

s x s y = U exp(i[4> x - 4> y \) 

= 3? exp(i[(p x - <f> y ]) exp(i[pirai(a;i - y±) + p 2 n 2 (x 2 - y 2 )\) , (33) 

where we have interpreted (j) x as the angle of the spin s x . 
The expectation value in the spin-wave limit becomes 

{SxSy)sW = 

Sm,n 2 W{n u n 2 ) {exp(i[(j) x - (j)y]))o fi cos(p 1 n 1 (x 1 - y x ) + p 2 n 2 (x 2 - y 2 )) 

E m ,n 2 W{n^n 2 ) 

(34) 

where (...)o,o denotes the expectation value in a system with vanishing bound- 
ary shift. Configurations with a winding (i.e. with a shift in 0) give only 
minor contributions; E.g. VF(1,0) = 3.487... x 10~ 6 for an L 2 lattice at 
P = 2/tt. 

We have computed (exp(i[<f) x — 4> y ]))o,o numerically, using the lattice prop- 
agator. To this end, we have used lattices up to L = 2048. For details of this 
calculation see the appendix. The results for < s x s y > were plugged into the 
definition (5) of the second moment correlation length. Extrapolating the 
finite lattice results to L — > oo gives 

& nd /L = 0.7506912... + 0.66737... {(5 - 2/tt) + (35) 
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Inserting ln l +c = n(/3 — 2/tt) for the critical trajectory, we obtain 

D 21 2430 

i 2nd /L = 0.7506912... + ^ * c ' + - ■ (36) 

Note that a similar result for the exponential correlation length on a 
lattice with strip geometry, i.e. an Lx oo lattice, can be found in the literature 
[23]: 

U P /L = 2(3 . (37) 
Inserting ln L 1 +c = ir((3 — 2/tt) into (37) gives 

( "' /L = \ + 1 tolTc + - (38) 

at the KT-transition. This prediction had been compared with Monte Carlo 
results in ref. [24] for lattice sizes up to L = 64. 
It is interesting to note that the limit 

\im J exp /L\ z=L/ , (39) 

where ^ e xp,oo is the exponential correlation length in the infinite volume limit 
in the high temperature phase, is exactly known for any z = L/^ exPtOD [25]. 
Note that this limit corresponds to the RG-trajectory that flows out of the 
point (x, z) = (0, 0), while the present study is concerned with the trajectory 
that flows into (x,z) = (0,0). 



4 Monte Carlo Simulations 

We have simulated the XY model at (3 = 1.1199, which is the estimate of 
ref. [7] for the inverse transition temperature and f3 = 1.12091 which is 
the estimate of Olsson [6] and consistent within error-bars with the result 
of Schultka and Manousakis [19]. For both values of f3, we have simulated 
square lattices up to a linear lattice size of L — 2048. The simulations 
were performed with the single cluster algorithm [26]. A measurement was 
performed after 10 single cluster updates. In units of these measurements, 
the integrated autocorrelation time of the magnetic susceptibility is less than 
one for all our simulations. 
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Table 1: Monte Carlo results for the helicity modulus T, the second moment 
correlation length over the lattice size ^2nd/L and the magnetic susceptibility 
X for two dimensional XY model on a square lattice of linear size L at (3 = 
1.1199. 



L 


T 


£,2nd/ L 


X 


16 


0.72536(7) 


0.79953(17) 


133.011(9) 


32 


0.70883(7) 


0.79231(18) 


452.114(31) 


64 


0.69785(7) 


0.78701(18) 


1536.58(11) 


128 


0.69001(7) 


0.78310(18) 


5220.99(36) 


256 


0.68400(7) 


0.77977(19) 


17729.9(1.2) 


512 


0.67926(6) 


0.77745(18) 


60185.8(4.0) 


1024 


0.67544(7) 


0.77532(19) 


204160.(15.) 


2048 


0.67246(10) 


0.77300(28) 


692146.(74.) 



For each lattice size and /3-value we have performed 5.000.000 measure- 
ments, except for L = 2048 were only 2.500.000 measurements were per- 
formed. We have used our own implementation of the G05CAF random 
number generator of the NAG-library. For each run, we have discarded at 
least 10000 measurements for equilibration. Note that this is more than what 
is usually considered as safe. On a PC with an Athlon XP 2000+ CPU the 
simulation of the L = 2048 lattice at one value of (3 took about 76 days. 

In table 1 we have summarised our results for the helicity modulus T, 
the second moment correlation length over the lattice size ^2nd/ L and the 
magnetic susceptibility x at /?= 1.1199. In table 2 we give analogous results 
at = 1.12091. 

First we fitted the helicity modulus T with the ansatz 

T = 0.63650817819 + const/ (In L + C) , (40) 

where const and C are the free parameters of the fit. Note that 0((lnL) 2 ) 
corrections that are due to e.g. the 0(z 2 ) contribution to T are effectively 
taken into account by the fit parameter C. Also corrections [3] proportional 
to In | lnL|/(lnL) 2 contribute to the value of C, since In | lnL| varies little for 
the values of L that enter into the fits. 

The results of the fits for f3 = 1.1199 are summarised in table 3 and for 
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Table 2: Same as table 1 but for (3 = 1.12091. 



L 


T 


{,2nd/ L 


X 


16 


0.72695(7) 


0.80044(18) 


133.174(10) 


32 


0.71059(7) 


0.79326(18) 


452.856(31) 


64 


0.69982(7) 


0.78888(18) 


1540.31(11) 


128 


0.69225(7) 


0.78464(18) 


5235.34(36) 


256 


0.68629(7) 


0.78157(19) 


17794.7(1.2) 


512 


0.68186(7) 


0.77951(19) 


60436.6(4.3) 


1024 


0.67826(7) 


0.77733(20) 


205185.(15.) 


2048 


0.67528(10) 


0.77547(28) 


696308.(75.) 



(3 = 1.12091 in table 4. For (3 = 1.1199 the x 2 /d-o.f. stays rather large even 
up to L min = 512. Also the value of C is increasing steadily with increasing 
L min . However this is not too surprising, since corrections that are not taken 
into account in our ansatz decrease slowly with increasing L. However, the 
results for const approach the theoretical prediction 0.318899454... as L m i n 
increases. For L m i n = 64 and 128, the x 2 /d-o.f. for f3 = 1.12091 is much 
larger than for (3 = 1.1199. However for L min = 256 it becomes about one 
for (3 = 1.12091. This should however be seen as a coincidence, since the 
value of const is increasing with L min and already for L min = 64 the value of 
const is larger than the value predicted by the theory. 

We conclude that our fit results are consistent with (3 = 1.1199 being the 
inverse transition temperature, while (3 = 1.12091 is clearly ruled out. One 
should notice however that fits with ansatze like eq. (40) are problematic, 
since corrections that are not included die out only very slowly as the lattice 
size is increased. 

Next we fitted the results for the second moment correlation length with 
an ansatz similar to that used for the helicity modulus 

t,2nd/ L = 0.7506912... + const / (In L + C) . (41) 

The results of these fits are summarised in table 5 for (3 = 1.1199 and ta- 
ble 6 for (3 = 1.12091. In contrast to the helicity modulus, we get a small 
X 2 /d.o.f. already for L min = 64. This might be partially due to the fact 
that the relative statistical accuracy of ^2nd/L is less than that of the helic- 
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Table 3: Fits of the helicity modulus at (3 = 1.1199 with the ansatz (40). 
Data with L = L min up to L = 2048 have been included into the fit. 



Lmin 


const 


C 


xVd.o.f. 


64 


0.2957(11) 


0.668(21) 


3.53 


128 


0.2988(17) 


0.740(37) 


2.67 


256 


0.3033(29) 


0.847(67) 


2.10 


512 


0.3097(52) 


1.01(13) 


1.77 


1024 


0.326(14) 


1.43(37) 





Table 4: Fits of the helicity modulus at (3 = 1.12091 with the ansatz (40). 
Data with L = L min up to L = 2048 have been included into the fit. 





const 


C 


X 2 /d.ol. 


64 


0.3382(13) 


1.201(14) 


16.56 


128 


0.3473(21) 


1.399(42) 


9.87 


256 


0.3616(36) 


1.724(79) 


1.03 


512 


0.3688(68) 


1.90(16) 


0.30 


1024 


0.377(16) 


2.09(40) 
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Table 5: Fits of the second moment correlation length of the lattice size 
£,2nd/L at (3 = 1.1199 with the ansatz (41). Data with L = L m i n up to 
L = 2048 have been included into the fit. 



Lmin 


const 


C 


X 2 /d.o.f. 


64 


0.2082(38) 


1.58(12) 


0.78 


128 


0.2086(58) 


1.59(20) 


1.03 


256 


0.2112(97) 


1.69(36) 


1.49 



Table 6: Fits of the second moment correlation length over the lattice size 
£,2nd/L at (3 = 1.12091 with the ansatz (41). Data with L = L min up to 
L = 2048 have been included into the fit. 



Lmin 


const 


C 


xVcl-o.f. 


64 


0.2435(47) 


2.26(14) 


2.24 


128 


0.2583(79) 


2.77(26) 


0.57 


256 


0.265(13) 


3.01(46) 


0.63 



ity modulus T. The result for const at (3 = 1.1199 is quite stable as L min 
is varied, and furthermore it is consistent with the theoretical prediction 
const = 0.212430... derived in this work. On the other hand, the fit results 
of const at (3 — 1.12091 are clearly larger than the theoretical prediction and 
furthermore the value of const is even increasing as L min is increased. These 
results are consistent with the analysis of the helicity modulus: While our re- 
sults are consistent with (3 = 1.1199 being the inverse transition temperature, 
(3 = 1.12091 is clearly ruled out. 

4.1 The magnetic susceptibility 

The magnetic susceptibility at the transition temperature is predicted to 
behave as 

X = const L 2 -" (lnL)~ 2r ... , (42) 
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with r = —1/16 and const depends on the particular model. This result can 
be obtained e.g. by integration of 

(s x s y ) oc iT 1/4 (In i?) 1 / 8 (43) 

given in ref. [3] for the correlation function, where R — \x — y\. Leading 
corrections to eq. (42) are due to the integration constant in eq. (18): 

X = const L 2 ~ v {lnL + C)- 2r ... . (44) 

In ref. [10] Irving and Kenna have simulated the same model as studied 
in this work on lattices up to L = 256. Using the ansatz (42), leaving r as 
free parameter, they find r = —0.023(10), which is about half of the value 
predicted by the theory. Later Janke [11] repeated this analysis for the XY 
model with the Villian action and lattices up to L = 512. He finds, also 
fitting with the ansatz (42), r = —0.0270(10), which is consistent with the 
result of Irving and Kenna. 

Here we shall check whether the value of r changes as larger lattice sizes 
are included into the fit. To this end, we only discuss the data for f3 = 1.1199. 
In table 7 we give results for fits with the ansatz (42), where we have taken 
— 2r as a free parameter. The x 2 /d.o.f. is very large up to L m i n = 256. For 
Lmin = 32 our results for — 2r is slightly larger than that of refs. [10, 11]. 
As we increase L min also — 2r increases. However, even for L min = 512, the 
result for — 2r is by more than 70 standard deviations smaller than the value 
predicted by the KT-theory. 

Next we checked whether this apparent discrepancy can be resolved by 
adding the leading correction predicted by the theory as free parameter to 
the fit. In table 8 we give our results for fits with the ansatz (44), where we 
have fixed — 2r = 1/8. We see that already for L min = 128 an acceptable 
X 2 /d.o.f. is reached. 

Finally we performed fits with the ansatz (44), where now also — 2r is 
used as free parameter. The results are summarised in table 9. The x 2 /d.o.f. 
becomes acceptable for L min starting from L min = 128. Now the fit results 
for — 2r for L min = 128 and 256 are consistent within the statistical errors 
with the theoretical prediction. 

We conclude that the apparent discrepancy with the KT-theory that was 
observed in refs. [10, 11] can be resolved by adding a correction term, which 
is predicted by the KT-theory, to eq. (42). 
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Table 7: Fits of the magnetic susceptibility at (3 = 1.1199 with the 
ansatz (42). Data with L = L min up to L = 2048 have been included into 
the Gt. 





const 


-2r 


X 2 /d.o.f. 


32 


0.9611(2) 


0.0699(1) 


382.5 


64 


0.9539(3) 


0.0741(2) 


119.2 


128 


0.9485(4) 


0.0772(2) 


35.7 


256 


0.9439(6) 


0.0798(3) 


5.2 


512 


0.9412(11) 


0.0812(6) 


1.5 



Table 8: Fits of the magnetic susceptibility at (3 = 1.1199 with the 
ansatz (44), fixing the exponent to the value — 2r = 1/8. Data with L = L m i n 
up to L = 2048 have been included into the fit. 



Lmin 


const 


C 


X 2 /d.o.f. 


8 


0.8121(1) 


4.423(9) 


307.2 


16 


0.8146(1) 


4.187(11) 


115.0 


32 


0.8170(2) 


3.953(14) 


32.5 


64 


0.8187(2) 


3.786(20) 


6.6 


128 


0.8197(3) 


3.690(28) 


1.5 


256 


0.8204(5) 


3.625(43) 


0.4 



Table 9: Fits of the magnetic susceptibility at (3 = 1.1199 with the 
ansatz (44). Data with L = L min up to L = 2048 have been included into 
the Gt. 



Lmin 


const 


C 


-It 


X 2 /d.o.f. 


32 


0.685(15) 


7.73(45) 


0.177(6) 


4.92 


64 


0.747(19) 


5.83(55) 


0.152(7) 


1.97 


128 


0.789(26) 


4.58(76) 


0.136(10) 


1.49 


256 


0.857(38) 


2.5(1.1) 


0.112(14) 


0.01 
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5 Summary and Conclusions 



We have studied the finite size behaviour of various quantities at the Kosterlitz- 
Thouless transition of the two-dimensional XY model. For the helicity mod- 
ulus T the value at the Kosterlitz-Thouless transition in the L — > oo limit 
and the leading logarithmic corrections to it are exactly known. Here, we 
have derived the analogous result (36) for the second moment correlation 
length over the lattice size £ 2 nd/ L'- 

6^ = 0.7506912... + ^^:+ ... . 

We have performed Monte Carlo simulations of the 2D XY model at f3 — 
1.1199 and j3 = 1.12091, which are the estimates of the transition tempera- 
ture of ref. [7] and ref. [6], respectively. Using the single cluster algorithm 
we simulated lattices of a size up to 2048 2 , which is by a factor of 5 2 larger 
than the lattices that had been studied in ref. [6]. Analysing our data for the 
helicity modulus T and the ratio ^2nd/ L we confirm (3 = 1.1199 as transition 
temperature, while (3 = 1.12091 is clearly ruled out. 

Fitting Monte Carlo data with the ansatze (40,41) is certainly a reason- 
able method to locate the transition temperature and to verify the Kosterlitz- 
Thouless nature of the transition. However one should note that the large 
values of x/d.o.f. of our fits and the running of the fit parameter C with 
the smallest lattice size L min that is included into the fits, indicate that sub- 
leading corrections that are not taken into account in the ansatze (40,41) are 
still large for the lattice sizes that we have studied. Since these corrections 
decay only logarithmically with the lattice size, it is difficult to estimate the 
systematic errors that are due to these corrections. 

Finally we studied the finite size scaling of the magnetic susceptibility. At 
the transition it should behave like x °c L 2 ~ v In L~ 2r with rj = 1/4 and r = 
— 1/16. However, fitting numerical data, the authors of refs. [10, 11] found 
r = —0.023(10) and r = —0.0270(10), respectively. Including larger lattices 
into the fits, our result for r moves toward the predicted value. Extending 
the ansatz to x °c L 2 ~ r ' (In L + C)~ 2r , where C is an additional free parameter 
consistent with the theory, the apparent contradiction is completely resolved: 
For a minimal lattice size L min = 256 that is included into the fit, we get 
r = -0.056(7). 
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7 Appendix: The correlation function at 

z = 

Here we compute the spin-spin correlation function for z — 0, i.e. for the spin 
wave approximation, for finite lattices with periodic boundary conditions. 

To this end let us first summarise a few basic formula on multi-dimensional 
Gaussian integrals as they can be found in text books on field theory. 

Our starting point is the generating functional 

| J D[0]exp^(0,A0)+^ =exp^(k,A- 1 k) ) j (45) 

where 

1 1 1 

2g(<f>, = 20 E AsyMy = 2g E [(</>* ~ ^+a) 2 + m2 €] (46) 

is the action of the Gaussian model on a square lattice and the partition 
function is given by 

Z = J D[4>]ew (—(<!>, A*)) (47) 



with 



/ m = n/#x . (48) 



For a square lattice with periodic boundary conditions A 1 can be easily 
obtained using a Fourier transformation: 

I e ip(x-y) 

l^i )xy T~9 / , ~~^> : n i 

L l p p z + m 2 
p 2 = 4 — 2 cospi — 2 cosj9 2 , (49) 
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where the pi, % — 1,2 are summed over the values {0, ...,L — 1} • (2n/L). 
Here we are interested in the massless limit m — > 0. Note that for Ex t = 
the contributions to (A;, A" 1 ^) from (pi,^) = (0,0) exactly cancel, while for 
J2 X kx 7^ 0, in the limit m — > 0, the right hand side of eq. (45) vanishes due 
to the divergent zero-momentum contributions to (k, A~ x k). Hence we get: 



lim n 4 / D[4>]exp 



-—(<!>, A<f>) + ik<f> 




liP(k,Ck) , if Ext — 
, otherwise . 



with 



G 



;r</ 



1 e ip ( x-1 ') - 1 



(51) 



Note that adding a constant to C xy does not change the result. Here we have 
chosen this constant such that C xx = 0. 

Now we are in the position to compute the two-point correlation func- 
tion (34) required for the computation of the second moment correlation 
length (5): 

(exp(i[(f) x - (p y ]))oo = exp [PC xy ] . (52) 

Due to translational invariance, it is sufficient to compute g{x) = C(o,o),x, 
for all lattice sites x. Employing the reflection symmetry of the lattice with 
respect to various axis the number of sites can be further reduced by a con- 
stant factor. Still, the direct implementation of eq. (51) would results in a 
computational effort oc V 2 for the calculation of £,2nd, where V is the number 
of lattice points. A more efficient method is discussed below. 

First we compute g(x) with x — (x±, 0) for x\ > 0: 



J/(*i>°) = £2 E Q(Pi) [e tpiXl - 1} 



(53) 



Pi 7^0 



with 



Q(pi) = E 



P2 



(54) 



I.e. these g(x) can be computed with an effort proportional to V. 
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Next we notice that g(x) satisfies Poisson's equation (see e.g. ref. [27] 
and refs. therein): 



Ag(x) - g(x - (1, 0)) - g(x + (1, 0)) - g(x - (0, 1)) - g(x + (0, 1)) = 
1 Ee ^ = f , if* = (0,0) 

L 2 ^ i -L , otherwise . v ' 

In principle, the remaining g(x) can now be computed recursively, using 
eq. (55). First one has to note that g(xi, 1) = g(xi, —1), where we identify 
L — 1 with —1, for symmetry reason. Therefore 

g( Xl , 1) = ^[Ag( Xl , 0) - g( Xl -1,0)- + 1,0) + L" 2 ] . (56) 
Then for x 2 > 1 one gets 

5-(xi,x 2 ) = 4g-(xi,x 2 - 1) 

-g{xi - 1, x 2 - 1) - g(xi + 1, x 2 - 1) - ^(xi, x 2 - 2) + Z+ 2 . (57) 

Unfortunately, rounding errors rapidly accumulate, and the recursion is use- 
less, at least when using double precision floating point numbers, for the 
lattice sizes we are aiming at. 

Instead, we have used an iterative solver to solve eq. (55). We imposed 
g(xi,0) = g(0,Xi) obtained from eq. (53) as Dirichlet boundary conditions. 
As solver we have used a successive overrelaxation (SOR) algorithm. With 
the optimal overrelaxation parameter, the computational effort is propor- 
tional to L 3 . We controlled the numerical accuracy of the solution by com- 
puting g(x) from eq. (51) for a few distances x. Since we could extract 
sufficiently accurate results for the limit L — > oo from lattice sizes up to 
L = 2048, we did not implement more advanced solvers like e.g. multigrid 
solvers. 
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